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Summary 

This preliminary investigation introduces the use of 
the Regier number as a flutter constraint criterion for 
aeroelastic structural optimization. Artificial Neural 
Network approximations are used to approximate the 
flutter criterion requirements as a function of the design 
Mach number and the parametric variables defining the 
aspect-ratio, center-of-gravity, taper ratio, mass ratio 
and pitch inertia of the wing. The presented approxi- 
mations are simple enough to be used in the preliminary 
design stage without a well defined structural model. 
An example problem for a low-speed, high-aspect-ratio, 
light- aircraft wing is presented. The example problem 
is analyzed for the flutter Mach number using doublet 
lattice aerodynamics and the PK solution method. The 
use of the Regier number constraint criterion to opti- 
mize the example problem for minimum structural mass 
while maintaining a constant flutter Mach number is 
demonstrated. 

Introduction 

The structural optimization of an airplane wing to 
satisfy flutter constraints is often an expensive process 
in terms of computer resource requirements. The typi- 
cal iteration of the design through the optimization pro- 
cess requires repeated determination of the unsteady 
aerodynamic forces and subsequent flutter solutions. 
While the solution of the flutter equations requires lit- 
tle computer resources, the calculation of the unsteady 
aerodynamic forces often requires significant computer 
storage and central processor time. One means to re- 
duce these costs, as discussed in references 1 and 2, is 
to use a simplified model of the unsteady aerodynamics. 
If the changes in the structure between adjacent itera- 
tions are small then the unsteady aerodynamic forces 
will be virtually unaffected and the aerodynamic model 
need not be updated at every iteration. A complete up- 
date might be made after, say, five iterations. Another 
simplified model, used in reference 3, assumes that the 
modal basis used to calculate the unsteady aerodynamic 
forces is unaffected during the optimization. This re- 
port presents an alternative approach to actually calcu- 
lating the unsteady aerodynamic forces and explicitly 
solving the flutter equations during the optimization. 
The method uses a flutter criterion to evaluate the flut- 
ter susceptibility of the wing. Certain properties of the 
wing are compared with the criterion to assess whether 
or not the wing has acceptable flutter characteristics. 
The use of a flutter criterion is particularly attractive 
during preliminary design when a variety of wings are 
under study, and high accuracy is not required. 

Although there are no flutter criteria that apply 
to all wings, there are criteria available that apply to 
many wings. For example, one flutter criterion is to 


have the sectional center-of-gravity forward of the sec- 
tional aerodynamic center-of-pressure. However when 
this criterion is used for the design of a wing, the re- 
sulting design has excessive weight. The criterion for 
the flutter constraint used in this report is based on the 
stiffness-altitude parameter, or Regier number, which 
depends on the stiffness of the wing and characteristics 
of the fluid in which the wing is operating, in partic- 
ular, the density and speed of sound. As the name 
stiffness- altitude implies, the value of the Regier num- 
ber increases as either stiffness or altitude is increased. 
This parameter has been used in reference 4 to corre- 
late the flutter results obtained from several hundred 
wind-tunnel, flutter-model tests. 

The proposed criteria is “If the Regier number for 
the wing being designed is greater than a reference 
value, the wing is flutter free.” The reference value is 
a calculated Regier number based on experimental flut- 
ter tests. If the Regier is substantially larger than the 
reference value, the wing has excess stiffness. This sug- 
gests that some material in the load bearing structure 
could be reduced ( i.e . reduce the stiffness) and thus 
reduce the structural mass. If the Regier number is less 
than the reference value, then the wing will flutter. An 
advantage of the proposed criterion over other criterion 
is that this proposed criterion is easy to calculate and 
simple enough so that it may be used early in the design 
process. 

The purpose of this paper is to discuss how a flutter 
constraint based on the Regier number can be used in 
optimizing an airplane wing. First the Regier number, 
including its background and meaning is introduced. 
This discussion is followed by the development of an 
artificial neural network (reference 5), or ANN, approx- 
imation for the reference Regier number data obtained 
from reference 4. Next, the procedure for using the 
Regier number in structural optimization is presented. 
Finally, the application of the Regier number criteria to 
an illustrative example is presented. 
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artificial neural network 

AR 

aspect ratio 
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speed of sound 
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reference length, wing semichord 
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structural damping coefficient 

height of aft or forward beam, Figure 8 

input of ANN scaling function, (Equa- 
tion 3) 

objective function 
reduced frequency 

ANN approximation for the AR correc- 
tion factor 

ANN approximation for the eg correc- 
tion factor 

ANN approximation for the A correction 
factor 

ANN approximation for the /i correction 
factor 

ANN approximation for the r a correc- 
tion factor 

semispan of wing 

Mach number 

sectional wing mass 

output of ANN scaling function, (Equa- 
tion 7) 

dummy argument of Equations 5 and 6 
Regier number 

Frueh’s modified Regier number 

Regier number approximation 

radius of gyration, normalized by 6 

Nonlinear transfer function used in ANN 

thickness of aft or forward beam wall, 
Figure 8 

vector of structural response variables 
velocity 

velocity index parameter 


W 


W a ,Wf 


Wi 


beam mass per unit length 

width of aft or forward beam, Figure 8 

weight parameter of neuron element 
(Equation 4) 

mass of the support structure 


x 


vector of structural optimization design 
variables 


Vj 

z 


to 

^ -irpb 2 

P 

Ps 

A 

A 

0, 


output of a neuron element (Equation 4) 

vertical displacement of the beam at the 
typical wing section 

frequency 

mass ratio 

mass density of fluid 

mass density of the support structure 

sweep angle 

taper ratio 

bias parameter of neuron element 
(Equation 4) 


subscripts: 

a pitch mode natural frequency 

a aft beam 

C conservative approximation 

d divergence condition 

£ best estimate approximation 

/ forward beam or flutter condition 

j j tfl neuron element 

max maximum value 

min minimum value 

r reference 

0 reference condition 

The Regier Number 

The Regier number, or stiffness-altitude parameter, 
is one of several nondimensional aeroelastic parameters 
that have been used to insure dynamic similarity be- 
tween models (reference 6). The stiffness-altitude pa- 
rameter was first suggested for use in displaying flutter 
data by Arthur A. Regier. Over time this parameter 
has become known as the Regier number, or 7Z. The 
expression for 71 is 


71 = 


Urby/jl 

a 


( 1 ) 


where co T is a reference frequency, b is the reference 
length, fi is the mass ratio, and a is the speed of 
sound. As is the case for the more familiar velocity 
index parameter 1 Vj, 71 can be derived by simplifying 

1 The velocity index parameter Vj is related to 71 by the 
following relationship: 


2 



the approximate empirical flutter formula expression 
given by Theodorsen and Garrick in reference 7. A 
modified version of TZ is presented by Frueh in reference 
8. Frueh ’s modified Regier number TZ is given by 


n = 


TZ 


( 2 ) 


where Cl is the lift curve slope. This expression would 
be attractive for use in correlating flutter data because 
TZ has a smoother variation with Mach number than the 
Regier number does. However, in many instances 
information is not available, especially from research 
flutter model tests, so it is often not possible to use 
TZ. It should be pointed out the Frueh’s expression 
can also be derived from Theodorsen and Garrick’s 
approximate empirical formula (reference 7) recognizing 
that C L a is 27 r in their formula. So, in effect, Regier’s 
formulation essentially assumes a constant lift curve 
slope. Recent research reports since about 1970 have 
not often used either TZ or Vj because these reports tend 
to be problem, or wing, specific. For a specific wing, 
presenting the stability boundary using TZ or Vj is not as 
useful as a plot of the stability boundary using dynamic 
pressure. However, in aeroelastic research reports prior 
to the 1970’s the focus of the research was in parametric 
studies and the presenting of the data using TZ or Vj is 
very common. 

1Z is also directly related to the parameters M and 
fj, that arise from a non-dimensional considerations of 
flutter. The term is the ratio of vibration velocity 
to fluid speed of sound and may be though of as a 
modified Mach number. The remaining term /x is 
the non dimensional ratio of mass of the body to the 
apparent mass of the fluid surrounding it. 

In reference 6 Regier describes 1Z as the square root 
of the ratio between the structural inertia force and the 
fluid force at Mach 1. Using the R,egier number to define 
the stability boundary stresses the importance of the 
ratio between the inertial forces and fluid forces as a 
primary factor in the governing equations of motion. 
Two dynamically similar models that have the same 
Regier number will have the same ratio of structural 
force to fluid forces at Mach 1. If the two models are 
also aerodynamically similar, then the flutter stability 
boundary will be the same when plotted in the R vs M 
plane. Figure 1 shows TZ plotted against Mach number 
for an typical wing. When constant dynamic pressure 
lines are plotted against M , they appear as radial lines 
through the origin as shown in Figure 1. The stable 
no flutter region is above the boundary; the unstable 
flutter region is below the boundary. Unlike Vj, TZ has 
the origin as an anchor point. The above properties 
makes TZ an attractive parameter for correlating flutter 
data from different wing configurations. 



Figure 1. Regier Number Stability Plane. 


For these reasons Harris (reference 4) chose TZ to cat- 
alog flutter data from a variety of wing configurations. 
Reference 4 is a summary of 341 experimental and the- 
oretical flutter studies presented as a series of flutter 
boundary plots of TZ vs M and correction factors, such 
as the correction factor for the aspect ratio, KaR vs AR. 
Figure 2 is a typical plot of 72- vs M taken from refer- 
ence 4 showing the flutter boundary and data points for 
low sweep conventional planform wings. 


s-r 

o 

rQ 


0 > 

'Eb 

<D 

« 


Si 

able 






Be 

m 

— r 
1 
1 
1 



6«0s?- 

!>is : 

Best 



Estimate 

Wk 

m 


■ 



Mach Number, M 

Figure 2. Regier Number vs Mach Number. 


Other data plots in reference 4 are used for conven- 
tional planform wing of moderate and high sweep, and 
delta wing designs. The reference 4 designated “conser- 
vative estimate” and “best estimate” curves are shown 
in Figure 2 as solid and dashed lines, respectively. The 
conservative estimate curve encompasses all the exper- 
imental data whereas the best estimate curve is similar 
to a mean fairing though the data. The data in Fig- 
ure 2 were adjusted by Harris to a nominal wing design 
for an average aspect ratio, center-of-gravity position, 
taper ratio, mass ratio, sweep angle, and radius of gyra- 
tion. In practice, the data from reference 4 is adjusted 
for the particular AR, eg, A, ^x, A and r a of the wing 
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design by using correction factors. The correction fac- 
tors are normally a function of a single parameter with 
the exception of the correction factor for p. This will be 
explained in the “The Regier Number Approximation” 
section of the paper. 

To determine if a wing design has an aeroelastic 
problem using reference 4, the 7Z of the design calcu- 
lated using Equation 1 is compared to the adjusted 7 Z 
value calculated by using the methods and data of ref- 
erence 4. A flutter condition is predicted if the 1Z calcu- 
lated using Equation 1 is less than the 7 Z generated by 
using reference 4. If the 1Z of the design is greater than 
the 1Z of reference 4, then addition weight savings can 
be achieved by removing some of the structural mass 
being used to generate “excess” this stiffness. 

The Artificial Neural Network 

Artificial neural networks (ANNs) are used in this 
report to approximate a number of curves in reference 
4. This section serves to give the uninformed reader 
a brief introduction to the ANN. Further information 
on ANNs can be found in reference 5. ANNs are not 
unique to the work presented in this report and several 
other functions could have been used to approximate 
the data. 


ANNs derive their name from trying to mimic the 
structure and functions of neural systems of living crea- 
tures. In general the network is composed of a number 
of signal processing elements, or artificial neurons. Fig- 
ure 3 shows the form of the ANN using the McCulloch- 
Pitts neuron model (reference 5). In the conventional 
use of ANNs, the data is approximated by a single large 
ANN with several inputs and outputs. Because the data 
in reference 4 have been separated into several single 
input/output functions the conventional ANN model is 
not used. All of the ANNs developed in this report have 
a single input and output variable. The input variables 
are initially scaled before being passed to each neuron 
elements as shown in Figure 3. The individual neuron 
elements reside in what is commonly called the “hidden 
layer.” The number of elements in the hidden layer de- 
termines the degree of accuracy of the approximation 
function. The more elements in the hidden layer the 
better the ANN is at approximating the original func- 
tion. The output of each of the hidden layer elements 
is passed to the “output layer” which is also a neuron 
element. The output of the output layer is again scaled 
before being used. 

The equations used for calculating the approxima- 
tion are: 

The input scaling function is 



I~M= ,y 'T 0 ' * 0.8 + 0.1 (3) 

\J-max *-min) 

where I sca i e d * s the output of the scaling function, and 
Imax and I-min are the input scaling parameters. 

The McCulloch-Pitts model for a single input neu- 
ron is described by 

Vj — *5*2 scaled 4" ®j) (^) 

where yj is the output, wj is the weight, and 0j is the 
bias of the neuron element. The function Si defines 
the type of limiting or nonlinear transfer function used 
in the neuron. The approximations generated in this 
report use the following two types of nonlinear transfer 
functions 

1 - e~ 2 P 
1 + e~ 2 P 

5 2 ( P ) = ^ (6) 

where p is the input to the function. The same type 
of nonlinear transfer function is used for all neuron 
elements of a particular ANN approximation. 

The output scaling function is 



0 = 


( Vj O-rnin ) 

R8 


*(o max d m tn) T 0 7 


(7) 
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where O is the output of the approximation, O max and 
O m i n are output scaling parameters, and yj is given by 
equation 4. The parameters that define the ANN are 
the minimum and maximum of the input and output 
scaling functions plus the weights and biases of the 
each of the neurons elements. Equations 3-7 are 
used to approximate each of the individual functions 
that are joined together to form the Regier number 
approximation as described in the next section. 

The Regier Number Approximation 

This section describes the approximations that are 
used to approximate the information contained in ref- 
erence 4. Reference 4 presents the design estimates of 
the flutter margins as a series of basic flutter boundary 
plots in the form of R vs M for wing configurations of 
specific values for AR, eg, A, y, A, and r a . Paramet- 
ric adjustment factors are presented as additional plots 
that are used to modify the basic boundary plots for 
variations in AR, eg, A, y, A, and r a • These adjust- 
ment factors are necessary to account for wing config- 
urations that are different from the base configuration. 
In Figure 4 the data and the ANN approximation func- 
tion are shown for the aspect-ratio adjustment factor 
Kar as a function of the wing aspect-ratio value AR. 
There are different basic flutter boundary plots for dif- 
ferent planform designs, such as high sweep or delta 
platforms. The approximations calculated in this re- 
port for the base flutter boundary are for subsonic, low 
sweep, conventional-planform wing configurations. 



Figure 4. ANN Approximation of I<ar- 

The Kar ANN approximation function is calculated 
using Equation 7. The Kar(AR) has a value of 1 
for AR = 2, which represents the aspect-ratio value 
of the base configuration wing. There are four other 
adjustment factors, K cg (cg ) for the sectional wing eg 


location, K\(\) for the taper ratio, I\ /I (y, A, M) for 
mass ratio, and K Ta (r a ) for the radius of gyration. 
Each of the adjustment factors are calculated using 
Equation 7. The mass ratio adjustment factor appears 
to be a function of y, A, and M but it is actually only 
a function of y. The parameters A and M determine 
which one of 6 different functions of y is used. The 
complete ANN approximation for the flutter boundaries 
are: 

v , = MM1 

C Ii AR (AR)K cg (cg)K x (\)K^, A, M)K r Jr a ) 

( 8 ) 

Tle(M) 

£ - K AR (AR)K cg (cg)K x (\)K^,A,M)K Ta (r„) 

( 9 ) 

where Rq an d Rg, each calculated using Equation 7. 
The functions Rq, equation 8, and R\, equation 9, are 
best and conservative estimate ANN approximations 
as a function of M, AR, eg, A, y, r a . The Rq 
approximation is based on the envelope of R whereas 
R *2 is based on the average value of R. 

Table 1 list the number elements in the hidden layer 
and the type of nonlinear transfer function (equations 
5 or 6) used for each function. The value of zero for 
the K Ta entry of Table 1 denotes no hidden layer. 
Table 2 lists the function coefficients for the ANN 
approximations. Note that the coefficients for K M entry 
of Table 1 is a single line in the table and is dependent 
on the value of A and M to determine which set of 
coefficients are used. As noted above and in Figure 
3 each ANN has input and output scaling functions. 
Table 3 lists the scaling coefficients used in equations 3 
and 7. 

Equations 8 and 9 describe simple approximations 
that can be used with a minimum of computational 
effort to develop a constraint function for use in a 
complex structural optimization procedure or to check 
a design for its flutter margin. 

Table 1. Regier ANN Configuration. 


ANN 

function 

Number of Elements 
in Hidden Layer 

Nonlinear Transfer 
Function 

r* e 

2 

Si 

R* c 

2 

Si 

Kar 

2 

s 2 

Keg 

2 

s 2 

K X 

2 

S 2 

Ky 

1 

s 2 

K Tq 

0 

s 2 
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Table 2. Regier ANN Coefficients. 


Function 

Hidden 

Wj 

Layer 

©> 

Output 

Wj 

Layer 

e? 

Kar 

-10.1802 

6.4287 

-2.8981 

-0.2088 


11.3170 

-1.6769 

2.5877 


Keg 

-8.8731 

4.6806 

1.8229 

-2.1408 


-12.3446 

0.9841 

5.6267 


K\ 

13.5425 

-1.5790 

-4.8732 

2.6204 


-9.4929 

4.8397 

1.7489 


AV 





M < .9 





A < 20° 

5.6802 

-2.1022 

-1.4161 

0.6581 

20° < A < 52° 

-6.1022 

1.4173 

2.7400 

-1.0061 

A > 52° 

6.0479 

-1.1682 

-3.4544 

2.3473 

M > .9 





A < 20° 

-6.2028 

1.0579 

2.7628 

-0.8023 

20° < A < 52° 

6.3106 

-0.8072 

-5.0643 

3.8784 

A > 52° 

5.3574 

0.6696 

7.5510 

-2.0054 

k Tq 



5.6931 

-2.8362 

n 

-1.3377 

-1.1461 

-0.3777 

0.6175 


1.4409 

-1.2542 

0.4905 


n 

1.3996 

-0.5984 

0.3697 

0.7787 


1.3784 

-1.0410 

0.1003 



Table 3. Regier ANN Scale factors. 


Input Output 


Function 

Max 

Min 

Max 

Min 

Kar 

2 

.2 

1.5000 

0.8993 

Keg 

60 

35 

1.7877 

0.8098 

K X 

1 

0 

2.2616 

0.9048 

Ap 

90 

10 

1.2390 

0.7512 

Kr a 

0.7 

0.3 

1.2630 

0.7321 

n* c 

1.8226 

0 

6 

-6 

n. 

2.6731 

0 

6 

-6 


Structural Optimization Using the 
Regier Number Constraint 

This section describes the use of the ANN approx- 
imation of the Regier number 1Z* in a structural opti- 
mization procedure for the design of flutter free wing 
while minimizing 

J = Wt (10) 

where J is the objective function and Wt is the mass 
of the support structure. Figure 5 shows the design 
optimization procedure used in this report. The proce- 
dure consists of two major independent computer pro- 


grams. These are a MSC/NASTRAN® 2 finite-element 
and sensitivity analysis program (reference 9) and, a 
second program denoted as Optimizer in Figure 5. 
The MSC/NASTRAN® analysis, denoted as (Sensi- 
tivity FEM Analysis), uses a finite-element model 
(FEM) to generate the coefficients of the sensitivity 
equations. The Optimizer program calculates the ob- 
jective function J , the gradient of the objective function 
and the constraint function g, for the optimization 
program CONMIN (reference 10). The CONMIN pro- 
gram is only used to solve the linear optimization prob- 
lem. The design variable equations that are generated 
by MSC/NASTRAN® have the form 

v - v 0 + ~( x - x q). ( 11 ) 

where v is a vector of structural design responses, vq is 
the value of the design responses at the start of the opti- 
mization step, is a matrix of sensitivity coefficients, 
xo is a vector of the value of the design variables at the 
start of the CONMIN optimization, and x is a vector of 
design variables used in CONMIN. The structural re- 
sponse vector v are functions that are computed within 
the MSC/NASTRAN® program. These responses can 
include weight, displacements due to specified load con- 
ditions, stress, and modal frequencies, etc. The el- 
ements of the design variable x are problem depen- 
dent and may include plate thickness, beam height, 
and flange width, etc. The matrix is computed 
by MSC/NASTRAN® using a finite-difference scheme. 
The sensitivities are transferred to the structural opti- 
mization program Optimizer as shown in Figure 5 and 
in more detail in Figure 6. Figure 6 shows the Opti- 
mizer’s internal flow as it cycles in the design space 
requiring the calculation of the objective function and 
constraint, or gradients of the objective function and 
constraint. As shown in Figure 6, the Optimizer pro- 
gram supplies the CONMIN program with the value 
of the objective function J, the constraint g, or the 
gradients of the objective function The CONMIN 
program calculates the gradient of g by using a finite- 
difference scheme. 

The Optimizer finds the value of the design vari- 
able x that minimizes the linear cost function of Equa- 
tion 10 subject to a non-linear inequality constraint 

iV-IKti ( 12 ) 

where 7 Z is Regier number calculated using Equation 1, 
and 1Z* is the required Regier number calculated from 


2 MSC/ is a registered trademark of The MacNeal-Schwendler 
Corporation. NASTRAN is a registered trademark of the Na- 
tional Aeronautics and Space Administration. 
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Figure 5. Design Optimization Iteration Procedure. 


either Equation 8 or 9. Hereafter Equation 12 will 
be referred to as the Regier design constraint. 7 Z is 
calculated using Equation 1 from the modal frequency 
(jj r , the wing mass m, constants dependent on the flight 
condition p and a, and if the shape optimization is 
being done, grid point locations of the FEM. The 
approximation 7Z* is calculated using either Equations 8 
or 9, and is a function of the independent variables AR, 
eg, A, A, /i, r a and M . With the exception of p and M , 
all of the independent variables for calculating 71 and 
71* must be calculated from the vectors v and x. The 
calculation of the Regier constraint from the vectors v 
and x is illustrated in the sample problem discussed 
in the next section. The values of p and M must 
be specified by the user. The CONMIN optimization 
program finds the value of x which will minimize the 
cost function, usually structural weight, subject to the 
Regier constraint defined by Equation 12. 

A complete design iteration is a single pass around 
the path shown in Figure 5, executing the Sensitiv- 
ity &: FEM Analysis and Optimizer programs once. 
At the end of the CONMIN optimization program, the 


design responses calculated with Equation 12 will not 
agree with the value calculated with the FEM analysis 
unless the two systems have converged. The design opti- 
mization procedure, Figure 5, must iterate between the 
Sensitivity &: FEM Analysis program and the Opti- 
mizer program until the design responses agree, or con- 
verge. A converged design is defined when the objective 
function and the design variables do not change dur- 
ing a design iteration. The reason for the disagreement 
between the CONMIN and MSC/NASTRAN® design 
responses is that the Sensitivity & FEM Analysis 
program is calculating the responses using a nonlinear 
FEM of the structure and the Optimizer program is 
using a model of the structure linearized about the ini- 
tial value of the design variables. At the end of the 
CONMIN optimization step the new values of the de- 
sign responses of the linear model may not agree with 
the responses of the nonlinear model at the optimized 
value of the design variables. This behavior can be seen 
in the results of the sample problem in the next section. 
If the design has not converged, a MSC/NASTRAN® 
analysis is done using a new FEM based on the new 
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value of x and the optimization of the linear structural 
model is repeated. The design optimization procedure 
is repeated until the design responses of the sensitiv- 
ity analysis agree with the value of the design response 
calculated with MSC/NASTRAN®, as shown in Figure 
5. 


CONMIN^ptimizatioiiProgran^J 



NO $ 



Figure 6. CONMIN Optimization Program. 


Example Problem Definition, Analysis, 
and Optimization Results 

This section describes an example problem for 
demonstrating the use of the Regier number constraint 
in the optimal design of a wing subject to a flutter 
requirement. The example was chosen to be complex 
enough to represent an actual aircraft wing. The fol- 
lowing model specifications were selected a priori: 

1. Shape optimization will not be done. 

2. The wing will have a high AR, AR > 5. 

3. The wing will be untapered, A = 1. 

4. The flow regime will be subsonic, and A = 0. The 
optimization will be for sea level flight condition. 

5. The support structure will be two simple beams. 

6. The dimensions and mass balance of the wing will 
be that of a typical light airplane wing that is in 
agreement with items 2-5 and designed to flutter 
at a subsonic M . 

As a result of the above specifications 7 Z* is a 
function of M, n, r a and eg. The optimization problem 
will be to design a minimum mass structure. The design 


Mach number is the flutter Mach number of the initial 
wing design. 
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Figure 7. Example Problem: Wing Planform. 


Structural Model 

The planform of the example wing is presented in 
Figure 7. The wing is rectangular in shape with a 
semispan and chord of 200 and 80 in., respectively. 
The structure of the wing is modeled as two uniform 
cantilever beams. The forward and aft beam locations 
are placed at 0.15c and 0.55c to represent a typical wing 
box structure. 

Each beam is modeled with 10 beam elements. The 
wing is constrained so that the aft beam is connect 
to the forward beam by rigid links. The forward 
beam is free to displace vertically and twist about its 
long axis. Each beam is carrying an amount of non- 
structural distributed mass that is proportional to the 
cross-sectional area from the the point midway between 
the beams to the leading or trailing edge of the wing. 
The forward beam has 0.16 lbm/in., and the aft beam 
has 0.44 lbm/in. of non-structural mass. To simulate 
the landing gear, concentrated masses of 16 and 23 lbm 
are located respectively on each beam at 20 and 40 in. 
span stations. To simulate fuel, concentrated masses of 
35 and 15 lbm are located respectively on each beam 
at 120 and 140 in. span locations. In Figure 7, these 
concentrated masses are represented by small circles on 
each beam. Initial mass and balance information for 
the example wing are given in Table 4. 

The design variables for the forward and aft beams 
are height hf and h a , width wj and w a , and wall 
thickness tf and t a of the rectangular beam section 
as shown in Figure 8. Table 5 summarizes the lower 
bound, initial, and upper bound values of the design 
variables. 


8 




Table 4. Example Problem: Wing Weights and Bal- 
ances. 


Element 

Structural 

Mass(lbm) 

Non-structural 

Mass(lbm) 

Total 

Mass(lbm) 

Forward Beam 

23.2 

32 

55.2 

Aft Beam 

23.2 

88 

111.2 

Fuel 


100 

100 

Landing Gear 


78 

78 

Total 

46.4 

298 

344.4 



Aft beam Forward beam 



Table 5. Example Problem: Wing Design Parameter 
Values. 


Parameter 

Lower Bound 
(in.) 

Initial Value 
(in.) 

Upper Bound 
(in.) 

Wf 

1 

2 

3 

h f 

2 

4 

5 

tf 

0.05 

0.1 

0.2 

W a 

1 

2 

3 

h a 

2 

4 

5 

ta 

0.05 

0.1 

0.2 


Flutter Analysis 

To determine the flutter Mach number design value, 
and to check the accuracy of the Regier constraint 


criterion, a separate flutter analysis is conducted at 
each optimization iteration. The flutter analysis is in 
no way a required step in the optimization procedure. 
The flutter Mach number My for the initial wing will 
define the M for the design calculations. This section 
gives the details of how the flutter speed is calculated 
using the MSC/NASTRAN® computer program. 

The first four vibration modes are used for the 
aeroelastic equations of motion. The vibration analysis 
results for these modes are presented in Table 6 and 
Figure 9. The first mode, 17 Hz, is the first bending 
mode. The second mode, 21 Hz, is the first torsion 
mode. The third mode, 46 Hz, is the second bending 
mode. The fourth mode, 58 Hz, is the second torsion 
mode. Both bending modes have significant amounts of 
torsion due to coupling between the plunge and pitch 
motion of the wing. 

Table 6. Example Problem: Wing Vibration Analysis. 


Mode 

Frequency 

(Hz) ' 

Generalized 

Mass 

Generalized 

Stiffness 

1 st bending 

17 

2.829 x 10 1 

3.271 x 10 & 

1 st torsion 

21 

1.769 x lO" 1 

3.090 x 10 3 

2 nd bending 

46 

1.350 x 10 1 

1.121 x 10 6 

2 nd torsion 

58 

6.760 x 10" 2 

8.980 x 10 3 


The unsteady aerodynamics are calculated using the 
doublet lattice method (reference 12). The aerodynam- 
ics planform, Figure 7, is divided into 10 spanwise and 
5 chord wise aerodynamic boxes for the double lattice 
procedure. Beam spline (reference 11) interpolation is 
used to calculate the aerodynamic grid deflection from 
the motion of the forward beam. Unsteady generalized 
aerodynamics forces, GAFs, are calculated for the val- 
ues of reduced frequency k = [.001, 0.1, 0.2, 0.3, 0.4, 
0.5, 0.6, 0.7, 1, 5] and M = [0 0.3 0.6]. Linear inter- 
polation is used to calculate the GAFs at intervening 
values of k and M. The PK solution method, refer- 
ence 11, is used to solve the dynamic equations of mo- 
tion at V = 600 to 12000 in. /sec in increments of 600 
in. /sec and M = [0, 0.2, 0.3, 0.4, 0.5, 0.6]. Figure 10 
presents a typical results of a PK analysis for M = 0.4 
for the above structural model with the design variables 
at their initial values. The data in Figure 10 indicates 
the I s * wing torsion mode going unstable at V = 5000 
in. /sec (M=0.37), c o = 19.5 Hz and I s * wing bending 
mode diverging at V = 9100 in. /sec, (M= 0.67). Note, 
that the predicated M for the data in Figure 10 is less 
that the specified M for the calculated GAF data. The 
true instability point, often referred to as the “match 
point”, is found by interpolating the data in Figure 10 
with other results from the PK solution method (M = 
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MODE 3, 46 Hz. MODE 4, 58 Hz. 


Figure 9. Example Problem: Vibration Modes. 

0.3). The interpolated match point, or stability condi- 
tion where p, M , and V are in agreement for a standard 
atmosphere, is Vj = 5030 in. /sec (Mj = 0.37), at uj = 
19.5 Hz and V d = 8590 in./sec (M d = 0.63). 

Regier Constraint Calculation 

The Regier constraint for use in the Optimizer, 
(Figures 5 and 6), is developed in this section. The 
elements of the vector v in Equation 11 for this example 
are the variables Wt, u> a , z a and Zj. The required 
parameters in Equations 8 or 9 for the calculation of 
the TZ* are M, p, r a , and eg. The purpose of this 
section is show how the these variables are calculated 
from the vector v. During the CONMIN optimization 
calculations, the design responses are approximated 
using Equation 11. For example, the deflection of the 
aft beam z a is approximated by 

c) z 

Za = Z a o + -q^( x ~ x o) (13) 

where z aQ is the value of the aft beam at the start of 
the CONMIN program. 

The Regier database in reference 4 is based on the 
sectional wing characteristics at the 75 percent span 


location. For this example problem the nearest grid 
locations occur at the 80 percent span location and 
these grid locations are used to calculate p, r a , and eg. 
The 80 percent span was used because for this problem 
the parameters p, r a , and eg do not change between 
the 75 percent and 80 percent span locations. 

The mass ratio for the example wing is 


Wt + 0.61 
rclb 2 p 


(14) 


where Wt is the mass of the support structure, l is the 
semispan of the wing, b is the wing semichord, and p is 
the density of the fluid. Note that this expression for p 
contains the mass of the beams and the nonstructural 
mass, but does not contain the mass of the landing gear 
or fuel. For the values of 7=200 in., 6=40 in. and p = 
1.39 x 10 -4 lbm/in. 3 Equation 14 becomes 


Wt + 120 
44.59 


(15) 


For the initial wing design p — 3.69. 
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The center-of-gravity eg is calculated by using 


^ x 100 
(16) 

where the sectional mass of the beams W a and Wj are 
calculated by using 

W x - p s (2t x (h x + w x ) - 4 (17) 

where p s is the mass density of the support structure, 
and the subscript x is either a or / to denote the aft or 
forward beam. For the initial wing design eg = 41.8. 

The radius of gyration r a parameter is calculated 
by using 

I (W a + 0.44) (0.1 - ea ) 2 + (Wy + 0.16) (0.7 - ea ) 2 
Ta = V Wa+W f + 0.6 

( 18 ) 

where ea is calculated by using the displacements zj 
and z a due to a torque applied to the wing tip. That 
is, 


0.05 ( Wa + 0.44) - 0.35 (Wy + 0.16) 
W a + W j + 0.6 



0 2000 4000 6000 BOOO 10000 


V in. /sec 



V in./sec 


Figure 10. Example Problem: Aeroelastic Analysis. 


ea = 0.1 — 0.8 — — — (19) 

z a + Zf 

For the initial wing design r a = 0.4. 

Figure 11 shows 1Z £ and 1Z £ as a function of M 
and the flutter analysis result (solid circle symbol) for 
the initial wing design. The data in Figure 11 indicate 
that the 1Z £ is a sufficiently accurate to be used in 
Equation 12 for 7 Z. That is, the approximation of “best 
estimate” stability boundary is essentially the same 
as the doublet lattice calculated stability boundary. 
Substituting Equations 1 for TZ and 9 for 1Z* into 
Equation 12, with Kar( 5) = 0.9029, Aa( 1) = 0.9028, 
and a = 13587 in. /sec yields the Regier constraint for 
an unswept wing with AR= 5 and A=1 as a function M , 
eg, p, and r a : 



Mach Number, M 


Figure 11. Example Problem: Regier Constraint Func- 
tion. 


1.2277 Z £ (M) u a ^/Ji q ( 20 ) Structural Optimization Results 

KMK cg (cg)K ra (r a ) 340 

The doublet lattice analysis predicts the initial ex- 
ample wing to flutter at M = 0.37. The objective of the 
where I<^( p) represents the adjustment factor ap- optimization is to design a minimum mass wing (be. 
proximation for M < .9 and A < 20°. minimize Wt) that will be flutter free at M — 0.37. For 
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M = 0.37 and 7Z £ — 0.621 Equation 20 becomes 

0 763 <fas/P_ Mil 

Kp(li)K cg {cg)K r „(r a ) 340 ' 

The design variable vector x in Equation 11 consists of 
the six variables w a , Wf, h a , hf,t a , and t f . CONMIN ’s 
solution to the problem is to set all of the design 
variables, except hf , to the lower design variable bounds 
and to find the value of hf that satisfies the constraint 
function defined by Equation 21. Figure 12 is the 
optimization history for the design variable hf and 
Table 7 lists the values of the design variables for the 
converged design. Figure 13 shows the optimization 
history of Wt. The final value of Wt is 16.6 lbm, a 
reduction of 29.8 lbm. 



Figure 12. Optimization Results: hf v s Iteration. 


Table 7. Example Problem: Final Wing Design Vari- 
ables Values. 


Design Variable 

Final Value 
(in.) 

Wf 

1 

h f 

4.474 

t f 

.05 

W a 

1 

ha 

2 

ta 

.05 


For the first and second optimization iterations in 
Figure 13 the discontinuous nature of the structural 
optimization results is quite noticeable. The first and 
second optimization iterations are bounded by the □ 
and A symbols. Convergence at the end of iteration 2 
between the linear optimization equations and the fi- 
nite element analysis equations is indicated in Figure 
13 by the symbol V, designating the third CONMIN 


optimization step, combining with the symbol A ) des- 
ignating the second CONMIN optimization step, to give 
the appearance of the symbol $. 



Figure 13. Optimization Results: Structure Mass vs 
Iteration. 


Figure 14 presents the history of the 7 Z and 7Z £ 
terms of the constraint function during the optimiza- 
tion. The solid line in Figure 14 is 71 and the dot- 
ted line is 7Z £ . Equation 20 is satisfied in Figure 14 
whenever the solid line is above the dotted line. At the 
beginning of the first iteration, denoted as iteration 0 
in Figure 14, the data indicates that the the structure 
has excess stiffness because 71 > 7 Z £ . During the first 
CONMIN optimization step, i.e. between iterations 0 
and 1 in Figures 12 - 14, the smallest hf is found that 
satisfies the nonlinear constraint described by Equation 
20. The reduction structural mass is indicated in Fig- 
ure 13. At the end of the CONMIN optimization step, 
the MSC/NASTRAN® analysis results are shown at it- 
eration 1 in Figures 12-14. The new constraint func- 
tion indicates that a further mass savings can be made 
by reducing the stiffness still further during the second 
CONMIN optimization step. The convergence of the 
optimization program in two iterations is confirmed by 
the third segment of the results being nearly flat in Fig- 
ures 12 - 14. 

In Figure 15 the flutter Mach number Mf predicted 
by the MSC/NASTRAN® flutter analysis is plotted 
against iteration number. The variation of Mf is 
considered to be within the accuracy of the calculation 
and the optimal design is considered to meet the design 
specifications. However, there is no guarantee that 
the optimal design will have the desired aeroelastic 
properties because the optimization only guarantees 
that 71 > 71*. 
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Figure 14. Optimization Results: 7Z and 71* vs Itera- 
tion. 



Figure 15. Optimization Results: MSC/NASTRAN® 

computed Mj vs Iteration. 

Only a single initial value case is done in this paper 
because the objective of the paper is the examination of 
the Regier number criterion and not the optimization 
problem. If the structural optimization problem is of 
interest, several different initial conditions would have 
to be tested determine the existence of other minimum 
of the objective function in the design space. 

Table 8 lists the execution times of the Sensitivity 
& FEM Analysis, Optimizer programs, and flut- 
ter analysis. The execution times listed in Table 8 are 
divided by the execution time of the flutter analysis 
program. Based on the data in Table 8, the Regier 
constraint is an efficient means for calculating an op- 
timal structure with flutter constraint. The execution 
time for the structural optimization using the Regier 
constraint will be less that the procedure described in 
reference 3 since the method in reference 3 requires a 
calculations of the GAF at each iteration of the opti- 
mization. 

Table 8. Example Problem: Program Execution Times. 


Program 

Relative Execution Time 

Sensitivity Analysis 

0.23 

Optimizer 

0.01 

Aeroelastic Analysis 

1.00 


Concluding Remarks 

A flutter constraint criterion based on the Regier 
number for aeroelastic structural optimization is intro- 
duced. The use of a constraint criterion makes it unnec- 
essary to calculate generalized unsteady aerodynamic 
forces between optimization cycles as do other methods 
used in aeroelastic design for flutter. 

The application of the method proceeds in the fol- 
lowing manner. Existing experimental flutter data in 
Regier number format are approximated by using Ar- 
tificial Neural Networks. This approximation is then 
used to develop the flutter constraint criterion. Next, 
the optimal flutter acceptable wing design is achieved 
by adjusting the design parameters to achieve the min- 
imum of the objective function, say structural weight, 
while satisfying the Regier constraint criterion. 

The procedure is illustrated by the application to 
an example problem. A simple, rectangular wing is op- 
timized for weight at a specified Mach number while 
achieving a desired minimum flutter condition. The 
data for developing the Regier number constraint are 
obtained from existing experimental results of similar 
configurations. The solution for the illustrative problem 
converged in two design cycles. Because of the simplic- 
ity of the Regier number approximation, little computer 
resources are required. 
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